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Abstract 

Rotationally coherent Lagrangian vortices are formed by tubes of deforming fluid elements 
that complete equal bulk material rotation relative to the mean rotation of the deforming fluid 
volume. We show that initial positions of such tubes coincide with tubular level surfaces of 
the Lagrangian-Averaged Vorticity Deviation (LAVD), the trajectory integral of the normed 
difference of the vorticity from its spatial mean. LAVD-based vortices are objective, i.e., remain 
unchanged under time-dependent rotations and translations of the coordinate frame. In the limit 
of vanishing Rossby numbers in geostrophic flows, cyclonic LAVD vortex centers are precisely the 
observed attractors for light particles. A similar result holds for heavy particles in anticyclonic 
LAVD vortices. We also establish a relationship between rotationally coherent Lagrangian 
vortices and their instantaneous Eulerian counterparts. The latter are formed by tubular surfaces 
of equal material rotation rate, objectively measured by the Instantaneous Vorticity Deviation 
(IVD). We illustrate the use of the LAVD and the IVD to detect rotationally coherent Lagrangian 
and Eulerian vortices objectively in several two- and three-dimensional flows. 


1 Introduction 

Coherent vortices still have no universal definition in fluid mechanics, but two main features of 
a possible definition have been emerging. First, vortices are broadly agreed to be concentrated 
regions of high vorticity. Some authors require this dominance of vorticity relative to other flow 
domains (McWilliams 1984, Hussein 1986), others expect it relative to the strain in the same 
domain (Okubo 1970, Hunt et al. 1988, Weiss 1991, Hua &: Klein 1998, Hua, McWilliams Sz Klein 
1998). Yet others compare vorticity to strain in the rate-of-strain eigenbasis (Tabor & Klapper 
1995; see also Lapeyre, Klein & Hua 1999, and Lapeyre, Hua & Legras 2001). 

Second, vortices are generally viewed as evolving domains with a high degree of material in¬ 
variance. Lugt (1979) writes that a vortex is a “multitude of material particles rotating around a 
common center”. McWilliams (1984) expects the vortex to “persist under passive advection by the 
large-scale flow”. Chong, Perry Sz Cantwell (1990) propose to capture vortices by finding spiraling 
particle motions in the frozen-time limit of the flow. Vortices are described as “highly impermeable 
to inward and outward particle fluxes” by Provenzale (1999), who requires small relative dispersion 
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within vortex cores (see also Cucitore Quadrio & Baron 1999). Chakraborty, Balachandar & Adrian 
(2005) proposes that both swirling motion and small particle separation should be distinguishing 
features of a vortex core. Haller (2005) views vortices as sets of trajectories with a persistent lack 
of Lagrangian hyperbolicity. Chelton et al. (2011) observe that nonlinear eddies (vortices with a 
rotation speed exceeding their translation speed) trap fluid in their interior and transport them 
along. Finally, in a similar geophysical setting, Mason, Pascual & McWilliams (2014) stress that 
vortices are “efficient carriers of mass and its physical, chemical, and biological properties”. 

The core of a coherent vortex is, therefore, broadly expected to be an impermeable material 
region marked by a high concentration of vorticity. What constitutes high vorticity is, however, 
subject to individual judgement, thresholding, and choice of the reference frame. It is therefore 
the material invariance of a vortex core that holds more promise as a first requirement in an 
unambiguous vortex definition. Indeed, the Lagrangian nature of a vortex can simply be assured 
by defining its boundary as a tubular (i.e., cylindrical, cup-shaped or toroidal) material surface. 

The challenging next step is then to select such a material surface in a way that it also encloses a 
region of concentrated vorticity. 

Unlike vorticity, materially defined vortex boundary surfaces are inherently frame-invariant, 
defined by a set of fluid trajectories rather than by coordinates or instantaneous scalar field values. 

In continuum mechanics terminology, a material vortex boundary must therefore be objective, i.e., 
invariant with respect to all Euclidean frame changes of the form 

X = Q{t)y + b{t), (1) 

where x G and y G denote coordinates in the original and in the transformed frame, re¬ 
spectively; Q{t) G 50(3) is an arbitrary rotation matrix; and b{t) G is an arbitrary translation 
vector (Truesdell & Noll 1965). Paradoxically, with the exception of the approach initiated by Ta¬ 
bor & Klapper (1995), none of the instantaneous Eulerian vortex criteria listed above are objective. 
Accordingly, they may only detect coherent structures after passage to an appropriately rotating 
or translating coordinate frame. Eor instance, the unsteady Navier-Stokes velocity field 

( xisin4t-|-X2(2-|-cos4f) \ 

xi(cos4t — 2) — X 2 sin4f j, (2) 

0 / 

is classified as a vortex by the Okubo-Weiss, Hua-Klein, Hua-McWilliams-Klein, and Chakraborty- 
Balachandar-Adrian criteria, as well as by the Q-criterion of Hunt et ah, the A-criterion of Chong, 
Perry and Cantwell, and the nonlinear eddy criterion of Chelton et al. In reality, (2) is a rotating 
saddle-point flow, with typical trajectories growing exponentially in norm. This instability, how¬ 
ever, only becomes detectable to these criteria after one passes to an appropriately chosen rotating 
frame (Haller 2005, 2015). Promisingly, if we simply impose the localized high-vorticity require¬ 
ment, the constant-vorticity flow (2) is immediately discounted as a vortex without further need 
for analysis. 

Selecting vortex boundaries as material surfaces ensures material invariance for the vortex, 
but any tubular material surface can a priori be considered for this purpose. Recent stretching- 
based variational approaches narrow down this consideration to exceptional material tubes that 
remain perfectly unfilamented under material advection (Haller &: Beron-Vera 2013, Blazevski & 
Haller 2004, Haller 2015). As an alternative, Farazmand &: Haller (2016) seek vortex boundaries 
as maximal material tubes along which material elements complete the same polar rotation over 
a finite time interval of interest. These approaches have proven effective in two-dimensional flows. 
They, however, rely on a precise computation of the invariants of the Cauchy-Green strain tensor 
along a Lagrangian grid, which requires the accurate numerical differentiation of trajectories with 
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respect to their initial positions. This presents a challenge in three-dimensional unsteady flows, for 
which the polar rotation approach additionally fails to be objective. Most importantly, however, 
Lagrangian strain-based approaches offer no link between material vortices and the expected high 
vorticity concentration, a defining feature of observed vortices. 

In summary, despite recent advances in vortex criteria and fluid trajectory stability analysis, 
a fully three-dimensional, computationally tractable and objective global vortex definition, with 
guaranteed material invariance and experimentally observable rotational coherence, has not yet 
emerged. Here we propose such a vortex definition and a corresponding vortex detection technique. 

Our approach is based on a recently obtained, unique decomposition of the deformation gradient 
into the product the two deformation gradients: one for a purely straining flow and one for a 
purely rotational flow (Haller 2016). This rotational deformation gradient, the dynamic rotation 
tensor, obeys the temporal superposition property of rigid body rotations, thereby eliminating a 
dynamical inconsistency of the classic polar rotation tensor used in classical continuum mechanics. 
The dynamic rotation tensor can further be factorized into a spatial mean-rotation component and 
a deviation from this rotation. The latter deviatory part yields an objective, intrinsic material 
rotational angle relative to the deforming fluid mass. 

We then define a rotationally coherent Lagrangian vortex as a nested set of material tubes, each 
exhibiting uniform intrinsic material rotation. Such a vortex turns out to be foliated by outward 
decreasing tubular level sets of the Lagrangian-averaged vorticity (LAVD). Additionally, we prove 
that the center of an LAVD-based vortex is always the observed attractor for nearby finite-size 
(inertial) particle motions in geostrophic flows. 

In the limit of zero advection time, our Lagrangian vortex definition turns into an objective 
Eulerian vortex definition; a set of tubular surfaces of equal intrinsic rotation rate. These surfaces 
are tubular level sets of the instantaneous vorticity deviation (IVD), providing a mathematical 
link between rotationally coherent Eulerian and Lagrangian vortices; the former are effectively 
derivatives of the latter. We illustrate our results on several examples, ranging from analytic 
velocity fields to time-dependent two- and three-dimensional models and observational data. 


2 Set-up 

We consider an unsteady velocity field v{x, t), defined on a possibly time-dependent spatial domain 
U{t) C over a finite time interval [toTi]- We assume that U{t) is invariant under the fluid flow 
generated by the velocity field (cf. eq. (8) below). Thus, U{t) is either a physical domain with 
an impermeable boundary, or U (t) is a material domain formed by a set of evolving trajectories of 
v{x, t). 

We write the velocity gradient Vu as 


'\Iv{x,t) = D{x,t)+ W{x,t), (3) 

where D = ^ (Vv -|- is the rate of stain tensor and IT = | (Vv — is the spin tensor. 

We recall that the vorticity w = V x u of the fluid satisfies 

1 

ITe = — 2 *^ ^ '^e e (4) 

We will also use the instantaneous spatial mean Co of the vorticity over U{t), defined as 


uj{t) 


/^(j)W(x,f) dV 
vol {U{t)) 


( 5 ) 
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where vol(-) denotes the volume for three-dimensional flows, and the area for two-dimensional 
flows. Accordingly, dV refers to the volume or area element, respectively, in U{t). 

Under general observer changes of the form (1), the spin tensor and the vorticity in the new y 
coordinate frame take the form 

W{y,t) = Q'^{t)W{x,t)Q{t) - Q^{t)Q{t), u;{y,t) = + q{t), (6) 

with the vector q defined uniquely by the relation QQ^e = \q x e for all e G (see, e.g., 
Truesdell &: Rajagopal 2009). Formula (6) shows that the spin tensor and the vorticity are not 
objective quantities: the eigenvalues and eigenvectors of W change in rotating frames, and so does 
the direction and the magnitude of oj. Thus, neither W nor a; is, by itself, suitable for defining 
distinguished material sets in the flow. This is because material sets are tied to evolving fluid 
particles without any reference to coordinates, and hence are inherently frame-invariant. More 
generally statements about the material response of a moving continuum cannot depend on the 
observer and hence should be objective (Gurtin 1982). 

Fluid particle trajectories generated by v{x,t) are solutions of the differential equation 

X = v{x, t), 


defining the flow map 


xo^ x{t;xo), fG[fo,U], (7) 

as the mapping from initial particle positions xq G U{to) to their later positions x{t;xo) G U{t). 
The assumed invariance of U{t) over the time interval [fojU] can now be conveniently expressed as 

Tl{U{t^)) = U{t) tG[to,h]. (8) 

We will also use the notion of a material surface, which is a smooth, codimension-one, time- 
dependent surface family M.{t) advected by the flow, i.e., M.{t) = Jf:^{M.{to)) t G [to)U]- 
The deformation gradient 

Fl{xo) = VFl{xo) (9) 

is a linear map, taking initial infinitesimal perturbations to the fluid trajectory x{t;xQ) at time 
to to their later positions at time t. Although often believed otherwise, is not objective: its 
eigenvectors and eigenvalues depend on the frame of reference (see, e.g., Liu 2004). Therefore, 
the invariants of F^^ do not provide an objective indication of the rotational component of the 
deformation. 

3 Finite material rotation from the Dynamic Polar Decomposition 

We seek to identify coherent Lagrangian vortices as the union of tubular material surfaces in 
which fluid elements exhibit the same bulk material rotation over a finite time interval of interest. 
Individual material fibers based at an initial point xq in a deforming continuum, however, all 
rotate around different axes and by different angles. In recent work (Farazmand Sz Haller 2016), 
we used the classic polar rotation angle (PRA) from finite strain theory to identify pointwise bulk 
material rotation in a moving fluid systematically. The use of the PRA, however, also leaves several 
challenges unaddressed, as we discuss in Appendix A. Most notable of these are an inconsistency of 
the PRA with experimentally observed dynamic rotation angles of spherical tracers in fluids, and 
its lack of objectivity in three dimensions. 

To address these challenges, we use here the recently developed Dynamic Polar Decomposition 
(DPD) to identify a dynamically consistent and fully frame-invariant rotational component in the 
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finite deformation of fluid elements (Haller 2016). This decomposition gives a unique, time-evolving 
factorization of the deformation gradient into the product of two deformation gradients: one for a 
purely rotational flow with zero rate of strain, and one for a purely straining flow with zero vorticity. 
Specifically, the unique right DPD of at xq can be written as 

Fl = OlMl, tG[to,h], ( 10 ) 

where the proper orthogonal dynamic rotation tensor = da^ait) is the deformation gradient of 
the purely rotational flow 

d = W {x{t; xo),t) a, (11) 

and the non-degenerate right dynamic stretch tensor Mj:^ = dbQb{t) is the deformation gradient of 
the purely straining flow 

h = Ol^D{x{f,xo),t)Olb. (12) 

The linear velocity field for a{t) in (11) is strainless because its coefficient matrix W is skew- 
symmetric. Similarly, the linear velocity field for b{t) in (12) is irrotational because its coefficient 
matrix 0^°D {x{t',XQ),t) Ol^ is symmetric. Unlike the classic polar rotation tensor (cf. Appendix 
A), the dynamic rotation tensor is dynamically consistent, i.e., satisfies the fundamental super¬ 
position property of solid-body rotations: 

0\^=OlOl^, s,t(^[to,ti]. (13) 

This follows because 0\^ is the fundamental matrix solution of a classical linear differential equation 
and hence satisfies the process property noted in (13) (cf. Dafermos 1971, Arnold 1978). In contrast, 
Mfg is the fundamental matrix solution of a non-classical linear differential equation with memory, 
i.e., with explicit dependence on the initial time to- Such fundamental solutions do not obey the 
process-property indicated in (13). The reason behind the dynamical inconsistency (35) of polar 
rotations is a similar memory effect in eq. (34). 

The decomposition in (10) is a right-type decomposition, i.e., the dynamic stretch tensor pre¬ 
cedes the dynamic rotation tensor from the right. Just as for the classic polar decomposition, a 
left-type version of the DPD is also available (Haller 2016). 

4 Lagrangian-averaged vorticity deviation (LAVD) 

Despite its dynamical consistency, the dynamic rotation tensor is not objective. Its frame- 
dependence is the consequence of the frame-dependence of the spin tensor W{x,t) appearing in 
the differential equation (11). The single remaining challenge out of those listed in Appendix A is, 
therefore, to identify an objective part of the rotation described by Ol^ which also preserves the 
dynamical consistency of Oj^. Below we recall further results from Haller (2016), and use them to 
address this challenge. 

The dynamic rotation tensor can further be factorized into two deformation gradients: one 
for a spatially uniformly rotating flow, and one for a flow that describes deviations from this uniform 
rotation. Specifically, we have 

= (14) 

where the proper orthogonal relative rotation tensor = daoOi{t) is dynamically consistent, 
serving as the deformation gradient of the relative rotation flow 

d = [W {x{t]Xo),t) — W (t)] a. (15) 
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Figure 1: The geometry of the intrinsic rotation angle ipl^^xo) along a material trajectory x(t) in a three- 
dimensional deforming continuum. 




In contrast, the proper orthogonal mean rotation tensor Q\^ = Di3^(3{t) is the deformation gradient 
of the mean-rotation flow 

$ = ^l°W{t)<!>ll3. (16) 

The mean rotation tensor is not dynamically consistent because (16) exhibits the same memory 
effect discussed for (12). 

The dynamic consistency of implies that the total angle swept by this tensor around its 
own axis of rotation is dynamically consistent. This angle called intrinsic rotation angle 

(see Fig. (1)), therefore satisfies 

V’Ua;o) = ifl{xo) + 'ipl^{xo), s,t £ [to,ti]- 

In addition, as shown in Haller (2016), '*/’*(,(xq) is objective both in two and three dimensions. 
In two dimensions, even the tensor itself turns out to be objective, not just its associated scalar 
field ifi^ixo). 

Using the results obtained in Haller (2016), the intrinsic dynamic rotation can be 

computed as 

<(xo) = ^LAVD^xo), (17) 

with the Lagrangian-Averaged Vorticity Deviation (LAVD) defined here as 

LAVD^i^(xo) := f |a;(x(s; xq), s) — a)(s)| ds. (18) 

Jto 

The objectivity of and LAVD can be confirmed directly from formula (6). Indeed, under a 
Euclidean observer change x = Q(t)y + b{t), the transformed vorticity uj{y,t) satisfies 

|w(y(s), s) - £5(s)| = |Q'^(s)w(x(s), s) + Q^{t)q{f) - {Q^{s)u{s) + {t)q{t)) \ 

= |Q'^(s) [a;(x(s),s) - a;(s)]| 

= |a;(x(s),s)-a;(s)| , (19) 

because the rotation matrix Q'^{s) preserves the length of vectors. We summarize the results of 
this section in a theorem. 
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Theorem 1. For an infinitesimal fluid volume starting from xq, the LAVD*^(xo) field is a dynami¬ 
cally consistent and objective measure of bulk material rotation relative to the spatial mean-rotation 
of the fluid volume U{t). Specifically, LAVD^^(xo) is twice the intrinsic dynamic rotation angle 
generated by the relative rotation tensor The latter tensor is obtained from the dynamically 
consistent decomposition 

Fl = (20) 

with the deformation gradient of a pure rigid-body rotation, and with the deformation gradient 
MtQ of a unique, purely straining flow. 

For detailed proofs of all statements in Theorem 1, we refer the reader to Haller (2016). Im¬ 
portantly, this theorem enables the extraction of an objective and dynamically consistent, material 
rotation component from the deformation gradient without carrying out the differentiation with 
respect to initial conditions in the definition (9) of K- 

5 Rotationally coherent Lagrangian vortices 

We now use the LAVD to identify objectively material tubes along which small fluid volumes 
experience the same bulk rotation over [to, ti] relative to the mean rigid-body rotation of the fluid. 
By Theorem 1, the time to positions of such material tubes are tubular level surfaces of the scalar 
function LAVDjJ(xo). By a tubular set, we mean here a convex, cylindrical, cup-shaped or toroidal 
set in three dimensions, and a closed convex curve in two dimensions. We require convexity for 
tubular surfaces, motivated by the near-circular cross section generally observed for stable vortices. 

If the gradient cIa;oLAVDjJ(xo) is nonzero along a tubular LAVD level surface, then this level 
surface is surrounded by a continuous, nested family of tubular level surfaces (Milnor 1963). The 
singular center of such a nested sequence of tubes, with inward increasing LAVD values, gives a 
definition of a Lagrangian vortex center. Similarly, the largest convex member of such a nested 
tube family defines the boundary of a Lagrangian vortex. We summarize these concepts in the 
following definition, with its geometry illustrated in Fig. 2. 

Definition 1. Over the finite time interval [to,H]-' 

(i) A rotationally coherent Lagrangian vortex is an evolving material domain V{t) such that 
V{to) is filled with a nested family of tubular level surfaces o/LAVDjJ(xo) with outward-decreasing 
LAVD values. 

(a) The boundary B{f) ofV{t) is a material surface such that B{to) is the outermost tubular 
level surface o/LAVDjJ(xo) in V{to). 

(in) The center C{t) ofV(t) is a material set C{t) such that C(to) is the innermost member 
(maximum) of the LAVDjJ(xo) level-surface family in V{to). 

We refer to the evolving positions C{t) of the tubular level sets 

^(^o) = {a^o £ U{to) ■ LAVD^J(xo) = const.} 

as a rotational Lagrangian Coherent Structure (rotational LCS), as indicated in Fig. 2. These LCSs 
give a foliation of the evolving Lagrangian vortex V{t) into tubes along which material elements 
complete the same intrinsic dynamic rotation fi\^{xo). 

Rotational LCSs, as well as Lagrangian vortices, their boundaries and centers are material ob¬ 
jects by definition. Therefore, their time t position is uniquely determined by Lagrangian advection 
of their initial positions: 

C{t) = (£(to)), B{t) = {B{to )), C{t) = (C(to)), t G [to,h]. 
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Rotationally coherent Lagrangian vortex 


Figure 2: Initial and current positions of a rotationally coherent Lagrangian vortexV(t); a Lagrangian vortex 
boundary B{t); and a Lagrangian vortex center C{t). Also shown within V{t) is a rotational LCS C(t), i.e., 
a material surface along which volume elements exhibit the same intrinsic dynamic rotation over the time 
interval As a conseguence, V{t) experiences only tangential material filamentation without global 

breakaway. 


Recent stretching-based definitions of Lagrangian vortices allow for no filamentation in boundary 
of the vortex (Haller &: Heron-Vera 2013, Blazevski & Haller 2014). In contrast, the LAVD-based 
definition of a rotational LCS allows for tangential material filamentation. The filamented part of 
the material surface, however, still rotates together with the LCS without global breakaway. 

The above definitions capture Lagrangian vortices with the simplest (i.e., convex) geometry at 
time to- More generally, one may allow for small tangential filamentation to be a priori present in 
the vortex boundary even at time to- This involves the relaxation of the convexity of £(to) and 
H(to) to material surfaces with small convexity deficiency, as discussed along with other numerical 
aspects in Section 9. 

In geophysical flows over a rotating planet, rotationally coherent Lagrangian vortices can di¬ 
rectly be computed from the flow induced in the curvilinear coordinate space instead of the curved 
surface of the planet (cf. Appendix B). By construction, the resulting vortices and their centers are 
invariant with respect to time-dependent rotations and translations within the space of curvilinear 
coordinates.^ With this approach, one simply computes the classic Euclidean vorticity of the lon¬ 
gitudinal and latitudinal coordinate speeds, as opposed to computing the vorticity in curvilinear 
coordinates. 

6 Rotationally coherent Eulerian vortices 

Over a short time interval [tQ,t + s] with |s| <C 1, we can Taylor expand the LAVD field as 

LAVD*+^(xo) = LAVD^xo) +IVD(x(t;xo),t) • s + O (s^) , (21) 

with the instantaneous vorticity deviation (IVD) defined as 

IVD(x, t) := |u;(x, t) — a)(f)| . (22) 

^Frame-invariance cannot be defined restricted to a curvilinear surface, as Euclidean frame changes take the 
observer off the surface. 
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Rotationally coherent Eulerian vortex 


Figure 3: A rotationally coherent Eulerian vortex V(t), with boundary B{t) and vortex center C{t). Also 
shown within V{t) is a rotational ECS E{t), i.e., a surface along which volume elements exhibit the same 
intrinsic dynamic rotation rate instantaneously at time t. No material coherence is guaranteed for the 
Eulerian vortex V (t) under passive advection. 


By the calculation (19), the IVD field is objective. By equation (21), IVD(x(t; xq), t) describes the 
rate of change of the LAVD field at an initial condition xq under increasing integration time. 

Using the IVD, we now introduce the instantaneous notion of a rotationally coherent Eulerian 
vortex by taking the limit to,ti —)• t in Definition 1. At a time t G [to,U]) such an Eulerian vortex 
is composed of tubular surfaces along which the intrinsic rotation rates of fluid elements are 
equal. Indeed, by formula (17), we have 'ipKx) = 2 lVD(x,t) = const, along these tubular surfaces. 
The following definition summarizes the details for this objective Eulerian vortex concept. 

Definition 2. At a time instance t G [tojU] : 

(i) A rotationally coherent Eulerian vortex is a set V{t) filled with a nested family of tubular 
level sets of IVD{x,t) with outwards non-increasing IVD values. 

(a) The boundary B{t) ofV{t) is the outermost level surface of IVD{x,t) in V{t). 

(Hi) The center C{t) ofV{t) is the innermost member (maximum) of the IVD(x,t) level-surface 
family in V (t). 

A rotational Eulerian Coherent Structure {rotational ECS) is then just a level surface 

E{t) = {x G U{t) : IVD(x,t) = Co(t)} (23) 

along which material elements experience the same dynamic rotation rate V'l(x) = Co{t)/2. Unlike 
rotational LCSs, individual rotational ECSs are instantaneous quantitates without a well-defined 
evolution, unless Co{t) is specifically selected as constant over time. We illustrate the geometry of 
Definition 2 in Fig. 3. 

Given that 

^LAVD^xo) = IVD(x(t;xo),t), 

rotationally eoherent Eulerian vortices are effeetively the derivatives of rotationally coherent La- 
grangian vortiees with respect to the length of the extraction time of the latter. Setting t = to in 
(21) and observing that LAVDj°(xo) = 0, we conclude that the B{to) initial positions of rotationally 


9 












Lagrangian vortex 
from LAVD under 
increasing extraction time 

Eulerian vortex 
from IVD at 
initial time 

Figure 4- ^ rotationally coherent Eulerian vortex is the limit of a rotationally coherent Lagrangian vortex 
under vanishing extraction time (i..e, integration) time. 



coherent Lagrangian vortex boundaries evolve precisely from their Eulerian coherent counterparts 
i3(to) as the Lagrangian extraction time ti — to increases from zero (see Fig. 4). 

Physically, these Eulerian vortices are built of tubular surfaces showing instantaneous coherence 
in the rate of their bulk material rotation. This instantaneous coherence in rotation rates generally 
does not imply sustained coherence in the finite rotation of material trajectories released from these 
surfaces. Furthermore, the interior of an evolving vortex V{t) is not a material domain. 

By the objectivity of the IVD, Definition 2 nevertheless gives an objective definition of an Eu¬ 
lerian vortex, its center and boundary. To our knowledge, no other objective, three-dimensional 
Eulerian vortex definition has been proposed in the literature. Specifically, none of the Eulerian 
criteria reviewed or proposed in Jeung & Hussein (1995), Haller (2005), and Chakraborty, Balachan- 
dar &: Adrian (2005) are invariant under time-dependent rotations and translations of the observer. 
Since truly unsteady flows have no distinguished frames of reference (Lugt 1979), non-objective 
Eulerian vortex criteria do not yield well-defined material vortices for unsteady flows. 

A relevant discussion can be found in Jeung & Hussein (1995) about the \uj\-criterion, by 
which |a;| must exceed a preselect threshold within a vortex. This approach is found intuitive but 
inadequate in Jeung &: Hussein (1995) for several reasons. We agree with this general assessment, 
because the |w|-criterion is threshold-dependent and not objective. In contrast, our rotationally 
coherent Eulerian vortex definition in Definition 2 is based on a threshold-independent and objective 
assessment of the level surface topology of Iw — a;| . 

Importantly, for flows with zero mean vorticity in their frame of definition, outermost convex 
tubular level sets of the vorticity magnitude |w| or of the enstrophy |ti;| coincide with rotation- 
ally coherent Eulerian vortices. Therefore, when properly interrogated, the vorticity and enstrophy 
distribution of zero-mean-vorticity flows do reveal objective structures that can be viewed as deriva¬ 
tives of rotationally coherent Lagrangian vortices. 

7 Rotationally coherent vortices in planar flows 

In a flow defined in the {xi,X 2 ) plane, the LAVD (18) takes the simple form 

t 

LAVD^xo) = y \uj 3 {x{s;xo),s) - u) 3 {s)\ ds, (24) 

*0 

with ids referring here to the X 3 component of the vorticity vector uj(x, t) = ( 0 , 0 , W 3 (x, t)), and with 
the mean vorticity 0 ) 3 ( 5 ) 

By Definition 1, a rotational LCS evolves over the time interval via advection from 

a closed and convex level curve of the LAVDjJ(xo) computed in (24). By the same definition. 
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Lagrangian vortex boundaries are outermost members of such rotational LCS families. Similarly, 
Lagrangian vortex centers are advected positions of isolated maxima of LAVD^|l^(xo). 

By Definition 2, a rotational ECS at time t is a closed and convex level curve of 

IVD(x,t) = \u:3{x,t) - a; 3 (t)| , 

around one of its local maxima. Accordingly, rotationally coherent Eulerian vortices are outermost 
members of such nested curve families with outwards non-increasing instantaneous IVD values. 
Eor two-dimensional flows only, Haller (2016) shows that the relative dynamic rotation angle 

1 

•= b / [‘^3(a;(s;a:o),s) -W3(s)] ds 
^ Jto 

is also an objective and dynamically consistent measure of rotation. It measures the net rotation 
angle generated by the relative rotation tensor around the X 3 axis, with sign changes in the rotation 
accounted for.^ 

While V’tQ(a^o)) as the total angle swept by the relative rotation tensor over the time interval 
is always positive, the sign of the angle is unrestricted. In vortical regions preserv¬ 

ing the pointwise the sign of the Lagrangian vorticity a; 3 (x(s; xq), s) — 0 ) 3 ( 5 ), contours of (/>^J(xo), 
contours of LAVD*)j(xo), and contours of the Lagrangian-averaged vorticity (LAV) 

fti 

■= u!3{x{s-,xo),s)ds 

Jto 

all coincide with each other, albeit generally correspond to different values of the underlying scalar 
fields.^ In such regions, the sign of (j)l({xQ) will also carry objective information about the direction 
of the relative rotation. 

Differences in the contours of LAVDJJ(xo), (j)l({xQ) and LAV^|)(xo) will arise in regions where 
the sign of a; 3 (x(s), s) — 0 ) 3 ( 5 ) crosses zero, i.e., near the boundaries of regions with a well-defined 
sign in their deviation from the mean vorticity of the flow. In such regions, connected level curves 
of LAVjJ(xo) can group together initial conditions with substantially different global rotational 
histories, as long as their final net rotations are equal. The use of LAVD^j) (xq) eliminates such 
coincidental agreement in the net rotation angles. 

8 Geostrophic Lagrangian vortex centers are attractors for iner¬ 
tial particles 

Consider a small spherical particle of radius tq and density Ppart in a geostrophic flow of density p 
and viscosity v. Under the /3-plane approximation, let / denote the Coriolis parameter (twice the 
local vertical component of the angular velocity of the earth). Applying a slow-manifold reduction 
to the Maxey-Riley equations (Maxey & Riley 1983) in the limit of small Rossby numbers, Beron- 
Vera et al. (2015) showed that the inertial particle motion satisfies 

x = v{x,t) + T{5-l)fJv{x,t) + 0{T‘^), 0 ^)’ 

where 

Ppart 9z/0 

^In contrast, 'iplg{xo) measures rotation about the evolving instantaneous axis of rotation of which points in 
the —X 3 direction when CJ 3 — 0)3 is negative. 

^Note that LAV// (xo) is not an objective quantity, but its level curves are objectively defined. 
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Figure 5: By Theorem 2, cyclonic Lagrangian eddy centers evolving from LAVD maxima act as observed 
attractors for small, positively buoyant particles, such as drifters and floating debris. 


Provenzale (1999) considered the Maxey-Riley equation in the same physical setting, but without 
a slow-manifold reduction. His second-order differential equation also included additional terms 
that either vanish along the /3-plane or appear at higher order in the reduced first-order equation 
(25) (cf. Beron-Vera et al. 2015 for more detail.) 

Remarkably, in the limit of vanishing Rossby numbers, cyclonic attractors for light particles 
{S > 1) and anticyclonic attractors for heavy particles (5 < 1) in eq. (25) turn out to be precisely 
the rotationally coherent vortex centers defined in Definition 1. The same Lagrangian vortex centers 
act as cyclonic repellers for heavy particles and anticyclonic repellers for light particles. We state 
these results in more detail as follows: 

Theorem 2. Assume that Xq G H(to) is the initial position of a rotationally coherent Lagrangian 
vortex center whose relative rotation keeps constant sign, i.e., 

sign[u) 3 {xl,t)- 03 {t)] = n{xl), te[to,ti], (26) 

for an appropriate sign constant ^ {“1 + !}• Then, for r > 0 small enough, the following 

hold: 

(i) In a cyclonic (/^(xq)/ > Oj rotationally coherent vortex, there exists a finite-time attractor 
(repeller) for light (heavy) particles that stays 0(t) close to the vortex center C{t). 

(a) In an anticyclonic (//(xq)/ < Oj rotationally coherent vortex, there exists a finite-time 
attractor (repeller) for heavy (light) particles that stays 0 {t) close to the vortex center C{t). 

Proof: See the Appendix C. 

Theorem 2 provides an independent, experimentally verifiable justification for defining vortex 
centers as in Definition 1 for geostrophic flows. Specifically, positively buoyant drifters or floating 
debris released well inside a cyclonic oceanic eddy will spiral onto the evolving Lagrangian vortex 
center identified from Definition 1 (Fig. 5). We will illustrate this effect using simulated inertial 
particle motion on satellite-based ocean velocities in Section 10.5. 

9 Numerical aspects 

Computing a rotation angle from the classic polar decomposition requires the computation of 
the deformation gradient K (cf. (33) in Appendix A). This is either achieved by the numerical 
differentiation of fluid trajectories with respect to their initial conditions, or by solving the equation 
of variations = Vn(x(t; to, xo),t)F(^, whose solutions typically grow exponentially. Either way, 
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computing polar rotation has the same long-time numerical sensitivity that arises in computing the 
invariants of the Cauchy-Green strain tensor. 

In contrast, computing rotationally coherent vortices based on Definition 1 only requires the 
integration of the normed vorticity deviation along fluid trajectories. This is the simplest to do 
simultaneously with trajectory integration, solving the extended system of differential equations 

= v{x,t), 

= |V X v{x,t) — V X v{t)\ , (27) 

over a grid of initial conditions over the time interval [toTi]- our experience, however, features 
of the LAVD^jj(xo) turn out to be sharper when the trajectory ODE x = v{x,t) is solved first, and 
the vorticity is subsequently integrated along trajectories. This is because adaptive ODE solvers 
make different decisions about time steps when the trajectory ODE is amended with the ODE for 
the LAVD held. This is especially so when both ODEs are solved simultaneously over large grids 
of initial conditions. 

The computational domain for solving (27) is just U{t) = U in case of a closed how with a 
hxed boundary. In open hows, the focus may be on vortices on a smaller domain. In that case, the 
domain should still be chosen large enough so that the averaged vorticity V x v{t) is representative 
of the overall mean rotation of the huid mass under study. In geophysical hows, this mean rotation 
is expected to be zero, which is conhrmed by our calculations even for domains of the size of a few 
degrees. 

For two-dimensional hows, we hrst identify local maxima of the LAVD held, then extract nearby 
closed LAVD level curves. In all our computations, we use the level set function of MATLAB for 
this purpose, and identify the closedness of a level curve by probing the output from this function. 
Dehnition 1 then requires the identihcation of the outermost convex LAVD level curve around an 
LAVD maximum as vortex boundary. This convexity requirement is conservative, ensuring that 
the material vortex starts out unhlamented at the initial time to- At later times, our approach 
allows for tangential hlamentation in the advected LCS due to local strain, but disallows large- 
scale hlamentation arising from differences in the bulk rotation along material hlaments. As a 
consequence, hlaments developed by rotational LCSs rotate together with the main body of the 
underlying material vortex. 

In actual computations, one reason to relax strict convexity for closed LAVD level surfaces is 
that they are numerically represented by discrete polygons. The more vertices such a polygon has, 
the more likely it is that the polygon is not convex, even if the approximated level curve is. A 
second reason for relaxing convexity is to remove the conceptual asymmetry of Fig. 2, allowing 
for small tangential hlamentation even in the initial positions of vortex boundaries. A third reason 
for relaxing strict convexity in multi-scale data sets is the presence of smaller-scale vortices near 
the perimeter of a larger-scale vortex. This necessitates the use of an appropriately coarse-grained 
notion for the boundary of the larger vortex. 

In two dimensions, addressing the hnite-grid, the initial tangential hlamentation, and the multi¬ 
scale challenges can be achieved by allowing for a small convexity dehciency in the LAVD level curves 
(Gonzalez & Woods 2008, Batchelor &: Whealan 2012). Here, we dehne the convexity deficiency 
of a closed curve in the plane as the ratio of the area difference between the curve and its convex 
hull to the area enclosed by the curve. In Fig. 6, we show cases of closed LAVD level sets with 
small convexity dehciency. In three dimensions, we require small convexity dehciency for the one¬ 
dimensional intersection of LAVD level surfaces with a family {Vi} of planes transverse to the 
expected vortex center curve C{to). 

In our computations, we set the convexity dehciency bound to 10“^ or lower, which, we hnd to 
produce robust results for the hrst two cases covered in Fig. 6. Applications to multi-scale data 


dx 

dt 
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LAVD level sets with small convexity deficiency 




Discrete approximation Level set has thin 
of a convex level set tangential filamentation 



Level set deformed by 
smaller-scale vortices 


Figure 6: Cases of closed LAVD level sets benefiting from the relaxation of convexity to small convexity 
deficiency. Shaded area indicates the area difference between the level set and its convex hull. (1) A convex 
LAVD level set with a non-convex approximating polygon arising from discretization. (2) Closed LAVD level 
set with minor initial filamentation in the tangential direction. (3) Closed LAVD level set that is only convex 
after a filtering of smaller-scale vortices. 


sets (last panel of Fig. 6) will likely require a somewhat higher bound. In general, increasing the 
convexity deficiency bound produces larger eddies that also tend to filament more. 

Small-scale local maxima of the LAVD function also arise due to numerical or observational noise 
in the velocity data. To eliminate the resulting artificial vortex candidates, we choose to ignore 
closed LAVD contours whose arclength falls below a minimal threshold. In a given application, 
this threshold should be chosen below the minimal vortex perimeter that is expected to be reliably 
resolved by the data set. 

We summarize the extraction algorithm for rotationally coherent Lagrangian vortices in two- and 
three-dimensional flows in the tables entitled Algorithm 1 and Algorithm 2 below. The extraction 
of their rotationally coherent Eulerian counterparts follows the same steps (2)-(4), but applied to 
the function IVD(x,f) instead of LAVD^|^(xo). 


Algorithm 1 Coherent Lagrangian vortex boundaries and vortex centers for 2-dimensional flows 
Input: A 2-dimensional, time resolved velocity field 

1. For a two-dimensional grid of initial conditions xq, compute the Lagrangian-Averaged Vor- 

ti 

ticity Deviation LAVDj|((xo) = /l^s (x(s; xq), s) — a; 3 (s)|(is. 

to 

2. Detect initial positions C{to) of vortex centers as local maxima of LAVD^|((xo). 

3. Seek initial vortex boundaries )3{to) as outermost, closed contours of LAVDjJ(xo) satisfying 
all the following: 

(a) I3(to) encircles a vortex center C{to). 

(b) S(to) has arclength exceeding a threshold /min- 

(c) i3(to) has convexity deficiency less than a bound dmax- 

Output: Initial positions of rotationally coherent Lagrangian vortex boundaries (closed curves) 
and vortex centers (isolated points) with respect to the time interval 
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Algorithm 2 Coherent Lagrangian vortex boundaries and vortex centers for 3-dimensional flows 
Input: A 3-dimensional, time-resolved velocity field 

1. For a three-dimensional grid of initial conditions xq, compute the Lagrangian-Averaged Vor- 

ti 

ticity Deviation LAVD^ilj(xo) = / |a; (x(s; xq), s) — a)(s)|(is. 

to 

2. Detect initial positions C{to) of vortex centers as local maximum curves (singular level sets) 
of LAVD^|((xo). This is best done by locating and connecting local maximum points Ci{to) of 
LAVDjJ(xo) in a family {Vi}^^ of N parallel planes. 

3. Seek initial vortex boundaries B{to) as outermost tubular level surfaces of LAVD^|((xo) whose 
intersections Bi{to) = B{to) n Vi with the planes Vi satisfy the following for i = 1,..., iV: 

(a) Bi{to) encircles a vortex center Ci(to). 

(b) Bi{to) has arclength exceeding a threshold /min- 

(c) Bi{to) has convexity deficiency less than a bound dmax- 

Output: Initial positions of rotationally coherent Lagrangian vortex boundaries (tubular surfaces) 
and vortex centers (curves) with respect to the time interval 


The procedure outlined in Algorithm 2 is mathematically well-defined by the level-surface topol¬ 
ogy of the LAVD. Step 2 of the algorithm describes one possible numerical extraction scheme for 
these level surfaces, using their intersections with planes transverse to the anticipated vortex center 
C{to). (A MATLAB implementation of this algorithm is available under https://github.com/LCSETH.) 
Complex flow geometries with multiple vortices and a priori unknown vortex orientations will re¬ 
quire more involved numerical approaches to level-set extraction. 

To implement Algorithms 1 and 2 in the forthcoming examples, we use a variable-order Adams- 
Bashforth-Moulton solver (ODE113 in MATLAB) for trajectory advection. The absolute and 
relative tolerances of the ODE solver are chosen as 10“® or higher. In Sections 10.3, 10.4, 10.5 and 
11.2, we use cubic and bilinear interpolation schemes for computing pointwise vorticity values for 
the IVD and LAVD functions, respectively. The lower-order interpolation suffices for the LAVD 
calculation, because trajectory advection smoothes out the LAVD contours in our experience. 

10 Two-dimensional examples 

10.1 Planar Euler flows 

On any solution of the two-dimensional Euler equation on the (xi,X 2 ) coordinate plane, the scalar 
vorticity value cj 3 (x,t) is preserved along trajectories. Therefore, the LAVD defined in (18) simpli- 
hes to 

LAVD4jj(xo) = (t-to) |w3(xo,to) -W3(to)| = (t-to)IVD(xo,to), 

Consequently, at time to, the boundaries and centers of all rotationally coherent Lagrangian vor¬ 
tices coincide with those of rotationally coherent Eulerian vortices. Specifically, B{to) = B{to) 
are outermost, closed and convex level curves of |a; 3 (xo,to) ~ ^ 3 (^ 0 )!) encircling a local maximum 
C{to) = C{to) of |w3(xo,to) - W3(to)|- 
If the planar Euler flow is steady, then 

= Vuj{x{t)) ■ v{x{t)) = 0. 
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Consequently, as long as the vorticity is not constant on open sets, vorticity contours coincide with 
streamlines. In that case, both the Lagrangian and the Eulerian vortices are bounded by outermost 
closed streamlines, with their centers marked by a center-type fixed point, as expected. 

Despite these formal LAVD calculations for planar Euler flows, one must remember: the eternal 
conservation of vorticity in these flows creates a high-degree of degeneracy for the rotation of fluid 
parcels. Indeed, despite the temporal and spatial complexity of a planar Euler velocity field u(x, t), 
a fluid element traveling along a fluid trajectory x{t) will keep its initial angular velocity (i..e, 
half of its vorticity at to ) for all times t. This property locks the initial conditions of initially 
rotationally coherent fluid parcels to the same LAVD level curve for all times, no matter how much 
these parcels separate or deform in the meantime. Specifically, even if material filaments break away 
transversely from a vortical region and undergo high stretching and global filamentation, they will 
eternally keep the exact same pointwise angular velocity they initially acquired near the vortex 
(see, e.g., the Eulerian vortex interaction simulation of Dritschel and Waugh 1992 for a striking 
example). 

This rotational degeneracy of the planar Euler equation is lost under the addition of the slight¬ 
est viscous dissipation, compressibility or three-dimensionality. Under any of these regularizing 
perturbations, LAVD-based vortex detection can be applied to reveal the Lagrangian signature of 
vortex interactions in the perturbed flow. We illustrate this in a viscous perturbation of the contour 
dynamics simulation of Dritschel and Waugh (1992) in Section 10.4. 

10.2 Irrotational vortices 

Although physically unrealizable, swirling flows with regions of zero vorticity are important the¬ 
oretical models. The classic irrotational vortex flow is a given by the two-dimensional, circularly 
symmetric velocity field 


vM = . (28) 

V J 

The vorticity of this flow is identically zero, implying LAVD^^(xo) = 0 and IVD{x) = 0 for any 
choice of to and t. This may seem to be at odds with the fact that all particles move on circular 
orbits, and hence line elements tangent to the trajectories exhibit one full rotation over the period 
of the circular orbit. 

However, the tangents are special line elements and are not representative of the overall bulk 
rotation of fluid elements. Material fibers transverse to the trajectory tangents all rotate at different 
speeds. The average rate of rotation for all line elements emanating from a given point is zero, as 
was already pointed out by Helmholtz (1858) in a related debate with Bertrand (1873) (cf. Truesdell 
& Rajagopal 2009, Haller 2016). Therefore, an infinitesimal fluid volume has no experimentally 
measurable bulk rigid-body rotation component in its deformation, as indeed demonstrated by 
Shapiro (1961). 

The irrotational vortex flow (28) is composed of a single, degenerate LAVD and IVD level set 
LAVD*^(xo) = IVD(x) = 0, as opposed to a nested family of codimension-one tubular level sets 
with outward decreasing LAVD and IVD values. Therefore, the flow (28) does not satisfy our 
Definitions 1 and 2, and hence does not qualify as a rotationally coherent Lagrangian or Eulerian 
vortex. A similar conclusion holds on the irrotational domains of the Rankine and the Lamb-Oseen 
vortices (Majda &: Bertozzi 2002), or more generally, in a recirculation region of any potential flow. 
We note that irrotational vortices are also either explicitly excluded by most systematic vortex 
criteria or fail the test for being a vortex according to these criteria. 
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10.3 Direct numerical simulation of two-dimensional turbulence 

We solve numerically the forced Navier-Stokes equation 

dtv + V • Vu = —Vp + V/S.V + /, V • u = 0, (29) 

for a two-dimensional velocity field v{x,t) with x = (xi,X 2 ) G U = [0,27r] x [0,27r]. We use a 
pseudo-spectral code with viscosity v = 10“^. We evolve a random-in-phase velocity field in the 
absence of forcing (/ = 0) until the flow is fully developed, then turn on a random-in-phase forcing 
(cf. Farazmand et al. 2013). We identify this latter time instance with the initial time to = 0, and 
run the simulation till the final time ti = 50. 

To construct the LAVD field, we advect trajectories from an initial grid of 1024 x 1024 points 
over the time interval [0, 50]. We integrated the vorticity deviation norm separately (as opposed to 
solving the combined ODE (27)), using 1200 vorticity values along each trajectory, equally spaced in 
time. Figure 7 shows the results superimposed on the contours of LAVDq°(xo). In this computation, 
we have set the minimum arc-length, Zmin = 0.3 and convexity deficiency bound dmax = 10“^. 

Figure 7 shows the Lagrangian vortex boundaries extracted from Definition 1 using Algorithm 
1 at the initial time to = 0. In Figure 7, we confirm the Lagrangian rotational coherence of these 
vortex boundaries by advecting them to the final time ti = 50. As expected, the vortex boundaries 
do not give in the general trend of exponential stretching and folding observed for generic material 
lines. Instead, they display only local (tangential) filamentation. The complete advection sequence 
over the time interval [0, 50] is illustrated in the online supplemental movie Ml. 

In Fig. 8 a, we show the time to positions of rotationally coherent Eulerian vortices in green, 
along with their Lagrangian counterparts in red. Passively advected positions of the Eulerian 
vortices are shown in Fig. 8 b, displaying substantial material fingering into their surroundings. 
While this lack of full material coherence over time is expected for the Eulerian vortices, some of 
them are impressively close to their Lagrangian counterparts at the initial time to, even though 
they are extracted just from an instantaneous analysis. At the same times, other Eulerian vortices 
without nearby Lagrangian counterparts disintegrate completely under advection, as seen in Fig. 
8 b. 

10.4 Interaction of two unequal vortices 

Motivated by the inviscid contour-dynamics simulation of Dritschel &: Waugh (1992), we seek to 
identify Lagrangian vortex cores in the interaction of two vortices of unequal strength (cf. our 
discussion in Section 10.1). We solve the Navier-Stokes equation (29) in two dimensions with 
viscosity v = 2 x 10“®. As initial velocity field, we let 

u(x,0) = i;i(x)-h U 2 (x), Vi{x)=Ti -— D /^i ) i = l,2, 

27rr^ ' 

where r? = jx —, and the matrix i? 7 r /2 refers counter-clockwise rotation by 7r/2. Each u* defines 
a vortex centered at a point Xi with a Gaussian profile. The parameters Tj and 5i control the 
strength and width of each vortex, respectively. In our simulation, we chose 

ri=0.01, (5i=0.05, xi = (7r,7r-F0.2), 

r2 = 0.12, (52 = 0.20, X2 = (7r,7r-0.2). 

The vortex described initially by vi is, therefore, weaker and smaller, and hence is expected to 
deform significantly under the strain field created by the second vortex initially described by V 2 - 
The computational domain is x G [0, 27r] x [0, 27r], discretized into 512 x 512 grid points. During the 
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(b) (c) 


Figure 1: (a) Three-dimensional contour plot 0 /LAVDg*^(a:o). Lagrangian vortex boundaries (outermost 
rotational LCSs) extracted from Algorithm 1 are shown in red. (b) Vortex boundaries shown in the flow 
domain at initial time t = 0. Shown in the background is the contour plot 0 /LAVDg*^(a;o) for reference, (c) 
Advected Lagrangian vortex boundaries at the final time ti = 50. (See the on-line supplemental movie Ml 
for the complete advection sequence of the vortex boundaries). 
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Figure 8: (a) Lagrangian (red) and Eulerian (green) rotationally coherent vortices at time to ^ the two- 
dimensional turbulence example, (b) Rotationally coherent Eulerian vortex boundaries and centers at time 
to, with the IVD(a;,to) field shown in the background, (c) The same objects advected passively to the final 
time ti = 50, with IVD(a;,ti) shown in the background. 
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Figure 9: (a) Rotationally coherent Lagrangian vortex boundaries (red) at time to = 0 in the vortex inter¬ 
action example. Also shown are two circles (black) approximating the Eulerian vortex boundaries inferred 
at time to from the vorticity distribution. The LAVD(^{xo) is shown in the background for reference, (b) 
advected position of the Lagrangian vortex boundaries and the two closed material curves at time t = 30. 
The advected LAVD field LAVDo^(x(2iQ]Q,xo)) is shown in the background. (See the on-line supplemental 
movie M2 for the complete advection sequence of the vortex boundaries). 


simulation, the vortices stay far enough from the boundary so that the boundary effects on their 
evolution are negligible. A relevant Reynolds number can be computed as Re = r 2 /iy = 6 x 10^ 
(cf. Kevlahan & Farge 1997). For this Reynolds number, the flow still remains fairly close to its 
inviscid limit, with the standard deviation of the vorticity from its initial value staying below 0.5% 
along fluid trajectories. 

We have run the computation from the non-dimensional time to = 0 up to t = 30. Over this 
time scale, we could preserve the smoothness of the velocity held without additional numerical 
effort. We applied Algorithm 1 with the tight convexity dehciency dmax = 10“^ and the arclength 
hlter Imin = 0.1. We show the rotationally coherent vortex boundaries obtained in this fashion in 
Fig. 9, along with two other closed material lines initiated around the vortex cores. The latter 
two circles approximate closely the Eulerian vorticity boundaries inferred from the instantaneous 
vorticity at time to = 0. They, therefore, mimic the the role of the initially circular vorticity jump 
contours advected in the inviscid simulation of Dritschel & Waugh (1992). 

As seen in Fig. 9, the stronger Eulerian vortex creates a sizable rotationally coherent Lagrangian 
vortex in the center. This vortex shows only tangential filamentation, as confirmed by its advected 
position at t = 30. The weaker Eulerian vortex has a substantially smaller coherent Lagrangian 
footprint that orbits around the larger vortex. 

This weaker satellite vortex gradually reaches a maximal elongated perimeter at about t = 24, 
then preserves its arclength in an approximate rigid-body rotation for the remaining 20% of the 
simulation time. No transverse filamentation occurs in this case either: the satellite vortex remains 
rotationally coherent in the sense of our Definition 1. 

Therefore, as Figure 9 illustrates, two unequal viscous vortices may interact strongly and still 
preserve their rotationally coherent material cores during a finite time interval of their interaction. 
Such coherent cores remain hidden in contour dynamics studies focused on the advection of Eulerian 
vortex boundaries inferred from large initial vorticity gradients. Indeed, the two black material 
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curves of Figure 9, mimicking the role of vortex-bounding inviscid vorticity contours, develop 
substantial transverse filamentation during the same time interval. 

10.5 Two-dimensional Agulhas eddies in satellite altimetry 

Here we illustrate the detection of rotationally coherent eddy boundaries in velocity data derived 
from satellite-observed sea-surface heights under the geostrophic approximation. In this approxi¬ 
mation, the satellite-measured sea-surface height r]{ip, 9, t) serves as a non-canonical Hamiltonian 
for surface velocities in {(p, 9) longitude-latitude {p, 9) coordinate system. 

The evolution of fluid particles satisfies 

*<*’■*•** = fl^/WcosO (31) 

where g is the constant of gravity, R is the mean radius of the Earth, and f{9) = 2Hsin0 is the 
Coriolis parameter, with H denoting the Earth’s mean angular velocity. The publicly available 
AVISO sea-surface height data base for ri{ip,9,t) is given at a spatial resolution of 1/4° and a 
temporal resolution of 7 days. 

We select the computational domain in the longitudinal range [—4°,9°] and the latitudinal 
range [—35°, —28°], which falls inside the region of the Agulhas leakage in the Southern Ocean. A 
Lagrangian analysis of coherent mesoscale eddies is particularly important in this context, as the 
amount of warm and salty water carried from the Indian Ocean to Atlantic Ocean has relevance 
for global circulation and climate (Beal et al. 2011). 

Recent two-dimensional Lagrangian studies of the Agulhas leakage used geodesic LCS theory 
to locate perfectly non-hlamenting (black-hole type) material eddies (Haller & Beron-Vera 2013, 
Beron-Vera et al. 2013). Here we use a relaxed notion of coherence, which allows for tangential 
filamentation, but not for global break-away of material from the eddy. The computational cost 
in the present method is substantially lower; the calculation of the deformation gradient and the 
search for limit cycles bounding the black-hole eddies are absent. 

We consider the AVISO data set ranging from the initial time to = November 11, 2006 to the 
final time H = to + 90 days. We select an initial grid of particles with stepsize Axq = 1/50°. As 
an additional filter for eliminating LAVD and IVD maxima due to resolution coarseness, we ignore 
from the start LAVD and IVD maxima that are closer than the submesoscale distance 0.2° (about 
20 km). As a complexity dehciency bound, we select dmax = 10“^. As an arclength threshold, we 
fix /mill = 27rr, with r = 20 km selected again as a lower bound on mesoscale structures reliably 
resolved by altimetry. 

Beyond executing Algorithm 1 to extract rotationally coherent in this setting, we also use this 
example to illustrate the predictions obtained from Theorem 2 for the attractor role of rotationally 
coherent Lagrangian vortex centers. Eigure 10a shows the rotationally coherent eddies, while 
Eigure 10b confirms that their boundaries only develop tangential filamentation under Lagrangian 
advection, as expected. This second plot also confirms that LAVD-based vortex centers are precisely 
the observed attractors for light particles released in cyclonic eddies, and for heavy particles released 
in anticyclonic eddies. 

We show a comparison of the initial positions of rotationally coherent Lagrangian and Eulerian 
vortices in Eig. 11a. Three of the Eulerian eddies (green) are close to Lagrangian eddies (red), and 
will accordingly show some coherence under advection by the end of the observational period in Eig. 
11c. The remaining Eulerian eddies show major filamentation and disintegrate under advection. 
We note the large number of false positives for coherent eddies based on the instantaneous Eulerian 
prediction at the initial time, even though this prediction is frame-invariant. 
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Figure 10: (a) Rotationally coherent Lagrangian vortices at time to = November 11, 2006, identified from 
Algorithm 1 using the contours of (xq) with T = 90 days. Shown in the background is the contour 

plot of LAVDI°'^'^(xq) for reference, (b) Initial (red) and final (black) positions of the Lagrangian vortex 
boundaries at time to + T, along with representative inertial particle trajectories. Heavy particles (blue) 
converge to the centers of anti-cyclonic (clockwise) eddies. Light particles (green) converge to the centers of 
cyclonic (clockwise) eddies. Here dueavy = 0.99, Sught = 1.01, rg = 1 m. (See the on-line supplemental movie 
M3 for the complete advection seguence of the vortex boundaries and the inertial particles). 
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Figure 11: (a) Lagrangian (red) and Eulerian (green) rotationally coherent vortices at time to in the satellite 
altimetry example, (b) Rotationally coherent Eulerian vortex boundaries and centers at time to, with the 
IVD(x,to) fie-ld shown in the background, (c) The same objects advected passively to the final time ti = to+T, 
with IVD(a:,ii) shown in the background. 
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In Appendix D, we also use this example to illustrate differences between LAVD-based vortex 
detection and two other objective Lagrangian tools for two-dimensional flows: the geodesic LCS 
approach of Haller & Beron-Vera (2013) and the ellipticity-time diagnostic of Haller (2001). 

We find that the geodesic LCS approach generally identifies similar material vortex regions 
as LAVD. By construction, however, geodesic vortex boundaries enclose perfectly non-filamenting 
vortex cores, and hence miss the rotationally coherent (but tangentially filamenting) outer annuli 
of Lagrangian vortices. The identification of the perfectly coherent, black-hole-type core by the 
geodesic LCS approach also comes with a higher computational cost and does not extend to three 
dimensions (cf. Appendix D for details). 

We also find that the ellipticity-time diagnostic highlights the general vicinity of the coherent 
material vortex regions, but offers no well-defined procedure for defining vortex boundaries. Focused 
on individual trajectory stability rather than global coherence, the ellipticity time has equally 
high values in some regions that do not remain coherent as a whole. Conversely, the ellipticity 
time indicates predominantly hyperbolic (i.e., non-vortical) behavior in the outer, tangentially 
filamenting parts of rotationally coherent vortices. 

In summary, despite their higher computational costs, neither globally uniform-stretching ma¬ 
terial lines nor elliptic fluid trajectories are able to provide the large, coherent material vortex 
boundaries obtained from the LAVD. 

11 Three-dimensional examples 

11.1 Strong Beltrami flows 

Strong Beltrami flows are steady flows whose vorticity field is a constant scalar multiple of the 
velocity field (Majda & Bertozzi 2002), i.e. 

Uj{x) = \v{x). 


In this case, formula (18) gives 


t 

LAVDiJxo) = |A| y |u(x(s)) - h| ds, v := j^v{x)dV. (32) 

Thus, for any strong Beltrami flow with A / 0, a rotationally coherent Lagrangian vortex bound¬ 
ary is a locally outermost, closed and convex level surface of |u(x(s)) — hi ds, the trajectory- 
averaged deviation of the velocity field from its spatial mean. For general 3D steady flows, level 
sets of asymptotically Lagrangian-averaged observables have been noted to approximate ergodic 
components (Budisic &: Mezic 2012). The new result here is that tubular levels sets of the finite¬ 
time average of the normed velocity deviation in strong Beltrami flows specifically define rotationally 
coherent Lagrangian structures (surfaces of constant net bulk rotation) in an objective fashion. 

By formula (32), rotationally coherent Eulerian vortices in strong Beltrami flows are composed 
of convex tubular level surfaces of the normed velocity deviation IVD(x) =|u(x) — v\. 

As an example, we consider the ABC flow whose velocity field is given by 

( A sin X‘s-\- C cos X 2 \ 

H sin xi-b A cos X 3 
C'sinx2-b Hcosxi / 

over the triply periodic domain Lf = [0, 27r]^. A direct calculation gives v{x) = 0 and A = 1. 
We select the parameter configuration A = 1, B = and C = For our Lagrangian 
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(a) (b) 

Figure 12: (a) Three representative rotationally coherent Lagrangian vortex boundaries in the ABC flow, 
obtained as outermost, convex tubular level sets of LAVDq® (xq). Black curves indicate trajectories launched 
on these boundaries, (b) Two representative IVD(a:o) level surfaces (with blue denoting the rotationally 
coherent Eulerian vortex boundary) in the single Eulerian coherent vortex obtained from the IVD field for 
the ABC flow for the parameter values considered here. 


advection, we select an initial grid of 200^ evenly distributed particles, which we advect over the 
interval [toj^i] = [0,50]. 

Figure 12a shows three representative Lagrangian vortex boundaries obtained as outermost, con¬ 
vex tubular level surfaces of LAVDq^ (a^o)- By the simplicity of this steady flow, relaxing convexity 
was not necessary (the convexity dehciency is zero). Also shown in Figure 12 are trajectory-segments 
launched along the Lagrangian vortex boundaries, illustrating their Lagrangian invariance. Figure 
12b shows the single IVD-based Eulerian vortex obtained for the ABC flow. While this Eulerian 
structure is near one of the Lagrangian vortices, the shape of the IVD-based vortex is not rep¬ 
resentative of the true Lagrangian vortex in this example. Furthermore, no IVD-based Eulerian 
vortices arise near the remaining Lagrangian vortices. This might seem puzzling first, but one must 
remember the spatial anisotropy of the ABC flow for unequal A, B and C parameter values. 

11.2 Three-dimensional Agulhas eddies in a data-assimilating circulation model 

Here we apply LAVD-based vortex extraction to a three-dimensional unsteady velocity held set 
obtained from the Southern Ocean State Estimation (SOSE) model (Mazloff, Heimbach & Wunsch 
2012). The domain of the data set is again in the area of the Agulhas leakage in the Southern 
Ocean, representing a three-dimensional extension of the Agulhas eddy extraction study of Section 
10.5. 

Our Lagrangian study covers a period of T = 120 days, ranging from to = 15 May 2006 to t = 
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12 September 2006. The selected computational domain in the South Atlantic ocean is bounded 
by longitudes [11°E, 16°E], latitudes [37°S,33°S], and depths [7,2000] meters. We compute the 
LAVD and IVD fields over a uniform grid of 350 x 350 x 350 points, and identify a rotationally 
coherent Lagrangian and Eulerian vortices using Algorithm 2. As in our two-dimensional turbulence 
example, we integrate the vorticity deviation norm separately (as opposed to solving the combined 
ODE (27)), using 1000 vorticity values along each trajectory, equally spaced in time. For the plane 
family Vi featured in Algorithm 2, we consider horizontal planes along nodes of the initial grid, 
starting from 28 meters below the sea surface level to eliminate noise due to boundary effects at 
the surface. In these planes we use the arclength threshold /min = 0.1 and the maximal convexity 
deficiency dmax = 10“^. The spatial mean vorticity a)(t), as the practically observed mean vorticity 
for a large enough fluid mass in the ocean, is taken to be zero. 

Figure 13a shows the initial position of a rotationally coherent Lagrangian eddy boundary 
(yellow) and its center (red), extracted as a level sets of LAVD^J(xo) by Algorithm 2. Also shown 
is a nearby LAVD level surface outside the eddy boundary, illustrating the complexity of the near¬ 
surface mixing region enclosing the eddy. Figure 13b gives a full view of the Lagrangian eddy, 
whereas Figure 13c shows the materially advected position of the eddy at the final time, 120 days 
later. As expected, there is mild tangential filamentation in the material eddy boundary, but 
strictly no break-away from the rotating water mass. Given the complexity of material mixing 
in the surrounding waters, this high degree of material coherence illustrates well the accuracy of 
LAVD-based vortex extraction. 

The corresponding rotationally coherent Eulerian eddy extracted at the same initial location, 
then materially advected for 120 days, is shown in Fig. 14. In the Eulerian computation, noise 
in the level surface computation is more moderate than in the Lagrangian case. As a result, the 
plane family Vi featured in Algorithm 2 can be selected to start from as high as 7 meters below 
the surface. 

This Eulerian eddy has about the same diameter near the sea surface as its Lagrangian coun¬ 
terpart, but maintains this diameter and reaches to substantially larger depths. Its vertical size 
is further increased under advection, with the bottom forming a sharp tip. Overall, the advected 
surface shows high ribbing and filamentation. Large-scale material break-away is absent in this 
example, but this cannot be guaranteed a priori for an Eulerian eddy (see our two-dimensional 
computations in Section 10.5). A comparison of Figures 13 and 14 suggests that the Lagrangian 
eddy forms a smooth, coherent center region that transports water without observable filamentation 
inside the vortical Eulerian feature. 

With a larger number of dedicated control points, the computation of LAVD level surfaces can 
be extended closer to the ocean surface. Figure 15 shows a higher-resolution computation of the 
initial and advected positions of the top slices of IVD-based and LAVD-based rotationally coherent 
vortex boundary surfaces. In this computation, the initial top slice of the LAVD-based vortex 
boundary is located only 15 m below the ocean surface, as opposed to the 28 m distance used in 
Fig. 13c. Even this Lagrangian boundary slice remains more coherent than the Eulerian one, but 
the Eulerian slice still performs well under material advection for this eddy and for this advection 
time. 

12 Conclusions 

We have given an objective (fully observer-independent) definition of a Lagrangian vortex as a 
set of material tubes in which fluid elements complete the same intrinsic dynamic rotation. This 
material rotation angle is obtained from the exact, dynamically consistent decomposition (20) of 
the deformation gradient. Remarkably, the intrinsic material rotation angle is expressible as the 
trajectory integral of the normed deviation of the vorticity from its spatial mean (LAVD). The 
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Figure 13: (a) Representative level surfaces of (xq) with to = May 15, 2006 and T = 120 days. 

The yellow surface is extracted using Algorithm 2, marking the vortex boundary for a mesoscale, rotationally 
coherent Lagrangian eddy, extending from 28 meters down to 646 meters in depth. The green surface is 
a nearby level surface of LAVD(‘’~'’^(a:o) outside the Lagrangian vortex region. The red curve marks the 
coherent vortex center, as defined in Definition 1. (b) Full view of the Lagrangian eddy boundary and its 
center at the initial time to (c) The advected eddy boundary and vortex center 120 days later, extending from 
49 meters down to 726 meters below the surface (See the on-line supplemental movie M4 for the complete 
advection sequence of the eddy boundary). 


intrinsic material rotation is, therefore, directly observable as the rotation of vorticity-meters in 
fluid experiments (Shapiro 1961), once the mean rotation reported by these devices is subtracted 
from the measurements. 

Locating a rotationally coherent material vortex does not require advection of high-density 
material grids, a generally taxing numerical procedure in Lagrangian coherence calculations. By 


27 







Latitude [°] 


Longitude [°] 

Latitude [°] 


Longitude [°' 


(a) 



(b) 



Figure 14: (a) The rotationally coherent Eulerian counterpart of the Lagrangian eddy shown in Figure 13a, 
ranging from 7 m down to 1,007 m below the see surface, (b) The materially advected position of the Eulerian 
eddy 120 days later, ranging from 15 m to 1,163 m below the sea surface (See the on-line supplemental movie 
M5 for the complete advection seguence of the eddy boundary). 


- LAVD top layer at t^^ 


Advected image of LAVD 
top layer at t^^+T 


- IVD top layer at t^^ 


Advected image of IVD 
top layer at t^+T 



Figure 15: Initial top slices (dashed line, at time to = May 15, 2006) and their advected positions (solid line, 
T = 120 days later) of the LAVD- and IVD-based rotationally coherent vortex boundaries. 
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construction, LAVD-based material vortices may show material filamentation, but even filamented 
material elements will rotate together with the vortex without global breakaway. The vortex inte¬ 
rior, therefore, shows no advective mixing with its environment, as we have illustrated on several 
several examples in two- and three-dimensional unsteady flows. 

Our approach also enables the extraction of vortex centers as singular level sets (local maxima) of 
the LAVD. We have proved that in two-dimensional geostrophic flows, these vortex centers coincide 
precisely with attractors of light particles in cyclonic eddies and for heavy particles in anticyclonic 
eddies. Indeed, we have found that numerically simulated light and heavy inertial particles show 
quick convergence to the appropriate moving vortex cores in a satellite-inferred, geostrophic ocean 
velocity field. On the same example, we have illustrated the advantages of LAVD-based coherent 
vortex detection compared to other objective Lagrangian vortex detection tools. 

The deviation of potential vorticity from its spatial mean is tempting to consider for a similar 
analysis, but such an approach would remain heuristic, driven purely by analogy. This is because 
potential vorticity is not objective and cannot be rigorously connected to intrinsic material rotation 
generated by the deformation gradient. 

Motivated by our results on LAVD, one might proceed by analogy and probe plots of La¬ 
grangian averages of arbitrary scalar fields for vortical features. As long as these scalar fields 
are non-degenerate and favorably initialized, material vortices should indeed have a footprint in 
the resulting plots due to the coherence of trajectories in the vortex interiors. For instance, the 
trajectory-averaged Okubo-Weiss parameter (Dresselhaus &: Tabor 1989, Prez-Muuzuri &: Huhn 
2013) and trajectory-averaged helicity (Prez-Muuzuri &: Huhn 2013) may also have lows and highs, 
respectively, near material vortices. Such extrema are, however, heuristic diagnostics, each signal¬ 
ing a different domain for a vortex, and each showing extrema in other flow regions as well. Their 
features are a consequence, rather than a well understood root cause, of material coherence (cf. 
Heron-Vera 2015). 

In contrast, we arrive here at LAVD-based vortices by solving an objectively posed coherence 
problem for material surfaces of equal bulk rotation. The LAVD then arises from an exact, dy¬ 
namically consistent decomposition of the deformation gradient into purely rotational and purely 
straining deformation gradients. The vortex boundaries and centers so obtained are sharply de¬ 
fined, do not develop global material filamentation, and remain invariant with respect to all possible 
Euclidean observer changes. 

We have also formulated an objective Eulerian definition of a rotationally coherent vortex: a do¬ 
main filled with tubular surfaces of constant intrinsic material rotation rate. These surfaces coincide 
with outward-decreasing tubular level sets of the instantaneous vorticity deviation (IVD). In some 
cases, we have found rotationally coherent Eulerian vortices and their centers to be surprisingly 
close to their Lagrangian counterparts. While such closeness does not hold in general, IVD-based 
vortices do provide a systematic and frame-invariant way to track coherent velocity features that 
are infinitesimally consistent in time with coherent material vortices. This makes these Eulerian 
vortices and vortex centers appropriate tools for a fully frame-invariant, automated vortex census 
in turbulent flow data. 
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13 Appendix A: Material rotation from the classical Polar De¬ 
composition 

Finite strain theory (Truesdell &: Noll 1965) infers a unique bulk material rotation for material 
volume elements between times to and ti from the polar decomposition 

= (33) 

of the deformation gradient. Here the rotation tensor is proper orthogonal, interpreted as the 
solid-body rotation component of the deformation. The right stretch tensor, , is symmetric and 
positive definite, obtained as the principal square root of the Cauchy-Green strain tensor C^. The 
right stretch tensor is interpreted as the stretching preceding the solid-body rotation represented 
by The tensors and Uj:^ are not objective, but the eigenvalues of are preserved under 
time-dependent rotations and translations (Truesdell &: Noll 1965). 

In a given frame of reference, RI^{xq) represents the unique rotation that gives the closest fit to 
the linear operator F^^(xo) in the Frobenius matrix norm (Gollub &: Van Loan 1983). In two and 
three dimensions, the action of RI^{xq) can be described by a single polar rotation angle (PRA). 
Well-defined up to multiples of 27r, PRA is the signed angle of rotation along the axis of rotation 
associated with the tensor R\^^{xo) 

In recent work, we used tubular and singular level surfaces of the PRA to visualize elliptic La- 
grangian regions and their centers, respectively (Farazmand &: Haller 2016). To our knowledge, this 
represents the first systematic approach to identifying Lagrangian coherence based on a synchrony 
in the net material rotation of infinitesimal volume elements. The level sets of PRA are objective in 
two dimensions, and have shown themselves to be accurate indicators of elliptic regions in unsteady 
flows with general time dependence. 

Nevertheless, several challenges remain that the classic polar decomposition, and hence also the 
PRA, are unable to address: 

1. PRA level sets are not frame-invariant in three dimensions. Therefore, PRA in three dimen¬ 
sions cannot be used to define elliptic LGSs, which are material (and hence fundamentally 
frame-invariant) surfaces. 

2. The computation of the PRA is based on the invariants of the Gauchy-Green strain tensor. 
This requires the accurate differentiation of the flow map (xq) with respect to initial 
conditions. This is numerically costly for extended times in large flow domains. 

3. No straightforward relationship exists between finite material rotation represented by the 
PRA and physical quantities (most notably, vorticity) used in Eulerian vortex identification. 
Indeed, only an involved link exists between R!^^ and the spin tensor W through the formula 

< = (»' - \k [R (R)"‘ - Kr' R] (Kf) K. (34) 

where dot refers to differentiation with respect to t (Truesdell & Rajagopal 2009). 

4. Between two fixed times to and ti, the polar rotation tensor R^^ represents the closest solid- 
body rotation to in the Frobenius norm. For intermediate times t G [to,G];the rotation 
family Rj^ provides no self-consistent solid-body rotation component for the evolving defor¬ 
mation gradients K . Indeed, one generally has 

-^<0 7 ^ s,t G (35) 
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which means that the rotation family does not satisfy the basic superposition principle 
of subsequent rigid-body rotations. Consequently, experimentally observed finite material 
rotation in fluids (visualized by small rigid-body tracers (Shapiro 1961) will differ from the 
PRA even in the simplest flows (see Haller 2016 for further discussion). 

Haller (2016) also discusses the rotational component yielded by Rl^ for two simple fluid-mechanical 
examples: irrotational vortices and parallel shear flows. In these examples, observed material 
rotation signaled by small inertial tracers (vorticity-meters) differs fundamentally from the rotation 
captured by Rj^. 

14 Appendix B: Rotationally coherent vortices in flow over a mov¬ 
ing surface 

Here we show how the LAVD and IVD can be computed over moving surfaces, such as the rotating 
Earth, without having to compute vorticity in the curvilinear coordinates. Consider first a three- 
dimensional spatial domain U(t) C that possibly also translates and rotates in time. We assume 
that U{t) admits a globally orthogonal parametrization 


ft:U C (36) 

a I—)• X, 

where a is non-dimensionalized. If the velocity field at points x GU{t) is denoted by v{x,t), then 

X = v{ff{a),t) = Dft{a)a + dtf{a), 

which induces the corresponding velocity field 

a = v(a, t) := [Dft{a)]~^ [v{ff{a),t) - dtf{a)] (37) 

in the parameter space U. The vorticity associated with the flow (37) in the parameter space is 
then given by 

oj{a,t) = Vq, X V, (38) 

with Vq, denoting the gradient operation in the orthogonal coordinates a. Theorem I is then 
applicable to the pull-back flow (37) in the parameter space with the vorticity 

As an example, we let U{t) denote a three-dimensional spherical shell region rotating with the 
Earth. A non-dimensional version of the classical spherical parametrization of the globe is given by 
a = (A/360°, V'/360°, p/i?), where A and denote longitudes and latitudes in degrees, p denotes 
altitude in kilometers, in and R denotes the radius of the Earth in kilometers. 

With the Earth modeled as a sphere placed of the x coordinate system and rotating with uniform 
angular velocity uq about the X 3 axis, we have the parametrization (36) in the form 

( 0:3 cos 27rQ;2 cos (27rai -|- vot) \ 

03 cos 27ra2 sin (27rQ;i -|- v^t) . (39) 

as sin 27rQ;2 / 

The velocity field v{a,t) in U{t) and the corresponding velocity field (37) in the space U of 
curvilinear coordinates are of the form 

/ vx{a,t) \ / vx{a,t)/[Ras cos {27rai + uot)] \ 

v{a,t) = vp{a,t) 1 , v{a,t) = vp{a,t)/ {Ras) I , (40) 

\ Vp{a,t) ) \ Vp{a,t) ) 
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respectively, Here {v\,v^,Vp) denote projections of the velocity v onto local unit vectors tangent 
to coordinate lines of the latitude, longitude and altitude. 

The parameter-space vorticity (38) can be computed from the velocity field v{a, t) in (40), then 
used in computing the IVD and the LAVD. For the special case of two-dimensional flows over a 
sphere, the p-component of v is zero, and p = R m all formulas above. 

15 Appendix C: Proof of Theorem 2 

The linearization of (25) along a particle motion Xp{t) starting from a position xq at time to gives 
the equation of variations 

^ = [Vv{xp{t),t) + T{d - l)fJVv{xpart{t),t) + C>(r^)] C- 

By Liouville’s theorem (Arnold 1978), the fundamental matrix solution Pl^{xo) of this linear system 
of ODEs satisfies the relationship 


detP* (^xq) = 


= ■l’to‘^3{Xpartis),s)ds + OiP) 


(41) 


By smooth dependence of the solutions of (25) on parameters, over a finite time interval and 
for small enough r, the inertial particle trajectory Xpart{s) is 0{t) C^-close to the fluid particle 
trajectory x{t]Xo) starting from the same initial position xq at time to. We thus have 

l^ 3 ixpart{s),s) = Ws(x(s; Xq), s) + 0{t), S € [fo,h]- 
Substituting this relation together with assumption (26) into (41), we obtain 

detPl (xq) = ^^ 2 ) 

with the two-dimensional LAVD field defined in (24), and the sign parameter^ defined in (26). 
Note that the planar form (24) of the LAVD applies in the present spherical flow setting because 
the /3-plane approximation is assumed in the derivation of the reduced Maxey-Riley equation (25). 

Over a finite time interval no unique classical attractor can be defined in a dynamical 

system. Indeed, by smooth dependence on initial conditions, any attracting trajectory has an open 
neighborhood filled with other attracting trajectories. What prevails from such an open set as a 
uniquely observed finite-time attractor is the trajectory that attracts the others at the strongest 
rate. 

Such a strongest attracting or repelling inertial particle motion starting from an initial po¬ 
sition xo is signaled by a local extremum of the function detP/j(xo) at Xq, i.e., by the relation 
VxQ [det Pj-^ (xq)] = 0. By (42), this extremum condition is equivalent to 

gM^(5-i)/LAVD(j(x5) _ i)/v^^LAVD*J(xS) + 0{t^)] = 0, 

which is in turn equivalent, for nonzero r, to an equation of the general form 

V.oLAVD*),(xS) + 0(t) = 0. (43) 

Assume now that Xq is a non-degenerate maximum point of LAVDjJ(xo), i.e., we have 

V,,LAVD*;(xS) = 0, det [v2^LAVD‘J(xS)] > 0. (44) 
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Then, by the implicit function theorem, for small enough r > 0, the equation (43) has a unique 
solution of the form 

xo(r) = Xq + 0(r). (45) 

Inertial particle trajectories starting from the initial position xo(t), therefore, prevail as the strongest 
finite-time attractors or repellers over the time interval These attracting and repelling tra¬ 

jectories remain 0(t) C^- close to fluid particle trajectories starting from the positions Xq. 

Now, on a large enough domain, the spatially averaged relative vorticity is approximately zero. 
(This already holds on the computational domain used in Section 10.5). Thus, the exponent in 
expression (42) can be written as 

r [fi{d - l)/LAVD*J(xS) + 0{t, /cug)] . 

When this expression is negative (positive) at the local maximum Xg of LAVD^^^ , then the fluid 
trajectory starting from Xg approximates the locally strongest finite-time attractor (repeller) of the 
inertial particle motion by (45). In other words, when ^{5 — l)f is negative (positive) at Xg, the 
Lagrangian fluid trajectory starting from Xg approximates the locally strongest attractor (repeller) 
over the time interval 

We conclude that in the limit of r —)• 0, cyclonic (/x/ > 0) Lagrangian eddy centers, as defined 
in Definition 1, are attractors for light (d > 1) particles and repellers for heavy particles (5 < 1). 
Likewise, for r —)■ 0, anticyclonic (pf < 0) Lagrangian eddy centers, as defined in Definition 1, are 
attractors for heavy (d < 1) particles and repellers for light (d > 1) particles. This completes the 
proof of Theorem 2. 


16 Appendix D: Comparison of LAVD-based vortex identification 
with other objective approaches 

16.1 Geodesic vortex detection 

Geodesic vortex detection seeks time to positions of Lagrangian vortex boundaries as outermost, 
closed stationary curves of the material-line-averaged tangential stretching functional 



A/(a^o('S),C‘o(®o(s))a;o(s)) 

— _^ 

V^(x'o(s),x()(s)) 


with C'|(j(xo) = [F/|j(xo)]^ T))j(xo) denoting the left Cauchy-Green strain tensor, and with r{s) 
referring to a parametrization of the closed curve 7 (cf. Haller & Beron-Vera 2013 and Haller 
2015 for details). On such stationary curves, the functional Q must necessarily have a vanishing 
variation; 

SQ{j) = 0. (46) 

This variational problem can be solved explicitly, with the solution depending on the eigenvalues 
Ai(xo) and eigenvectors ^i(xo) of Cfg(xo), defined and indexed as 

l^il = 1) ^ = 1)2; 0 < Ai < A 2 , Cl -L C 2 - 


Using these quantities, all closed curves solving (46) can be expressed as limit cycles of the 
autonomous differential equation family 


x'o = 


'A 2 -A 2 

A 2 — Ai 




A 2 — Ai 




(47) 
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for some value of A > 0 and for some choice of the sign in ±. This constant A turns out to be precisely 
the factor by which any subset of a trajectory of (47) will be stretched under the flow map 
Outermost members of nested limit cycle families of (47) are, therefore, locally the maximal closed 
curves in the flow that stretch uniformly (i.e., without filamentation). The geodesic theory of elliptic 
LCSs developed by Haller & Beron-Vera (2013) defines coherent Lagrangian vortex boundaries to 
be these outermost limit cycles. These boundaries are objective by the objectivity of the invariants 
of C'ljj(xo). An automated detection algorithm for geodesic vortex boundaries is given by Karrasch 
et al. (2014). Geodesic vortex detection has no direct extension to three-dimensional flows, but a 
variational approach for nearly-uniformly-stretching material surfaces is available (Ottinger et al. 
2015). 

16.2 Ellipticity—time diagnostic 

Haller (2001) studies the finite-time stability of a fluid trajectory x{t; to, xq) in a frame aligned with 
to the eigenvectors of the rate-of-strain tensor D{x, t) along the trajectory. A topological argument 
shows that the trajectory has an instantaneous elliptic stability type if throughout the time interval 
of interest, either D vanishes or the strain-acceleration tensor 

M = b + 2DVv, (48) 

is indefinite on the zero rate of strain set Z = |a G ; {a, Da) = O}. (In (48), dot refers to the 
material derivative.) The ellipticity time Te{ti,to,xo) for a trajectory released from xq at time to 
is then defined as the percentage of time within [to, ti] over which the trajectory has instantaneous 
elliptic stability. 

The scalar field Te{ti,to, xo) is an objective, pointwise indicator of fluid trajectory stability. It 
does not offer a strict dehnition of a material vortex boundary, but its high values indicate the 
general location of material vortices. As shown in Haller (2001), Te{ti,to,xo) can equivalently be 
defined as the percentage of time over which the trajectory x(t; to, xo) is elliptic in the sense of the 
Okubo-Weiss criterion, applied in a frame co-rotating rate-of-strain eigenbasis. A similar ellipticity 
time diagnostic can be dehned in three dimensions (Haller 2005), but this extension is no longer 
related to other known instantaneous vortex criteria. 

To be experimentally verihable, a vortex criterion based on coherent rotation of fluid elements 
should have a direct relation to observable mean material rotation visualized by small inertial tracers 
with an attached arrow (vorticity meters). Experiments show that this mean material rotation has 
an angular velocity that is precisely one half of the local vorticity (Shapiro 1967). In contrast, the 
relative vorticity observed in strain basis theoretical remains invisible under all possible Euclidean 
observer changes from the lab frame. Indeed, the coordinate change to the pointwise differing 
rate-of-strain eigenbases would require a spatially nonlinear rotation tensor Q{x,t) in (1), under 
which (1) no longer describes a physically meaningful observer change. 

16.3 Comparison with LAVD on the Agulhas leakage data set 

Both the geodesic and the ellipticity-time approach require more computational effort than the 
LAVD approach. Eor the geodesic method, limit-cycle families of a vector field composed of the 
invariants of the Cauchy-Green strain tensor must be computed with high accuracy, requiring the 
accurate differentiation of trajectories with respect to their initial conditions. Eor the ellipticity-time 
approach, the time-derivative of the rate-of-strain eigenbasis must be determined along trajectories 
wit high accuracy. 

Both approaches are also more stringent than the LAVD approach, requiring uniform stretch¬ 
ing (geodesic method) or a lack of trajectory-level instability (ellipticity-time diagnostic) along 
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Figure 16: (a) Ellipticity time field Te{ti,to,xo) for the Agulhas leakage data set of Section 10.6 with to = 
November 11, 2006 and ti = to + 90 days. Rotationally coherent Lagrangian vortex boundaries (red) are 
superimposed for reference, (b) Geodesic vortex boundaries at time to for the same data set, with the LAVD- 
based vortex boundaries (red) superimposed. 


the Lagrangian vortex boundaries. By their objectivity, both approaches should generally capture 
the same Lagrangian vortex region as the LAVD approach, but are expected to yield tighter vor¬ 
tex boundaries because of their more stringent definitions of coherence and trajectory stability, 
respectively. 

Figure 16 shows the results from these two alternative approaches for the data set analyzed in 
section 10.5, with the LAVD-based vortex boundaries superimposed in red. For the purposes of this 
comparison, we have used the numerical implementation of the geodesic eddy detection method 
described by Hadjighasem &: Haller (2016). 

As seen in Fig. 16a, the pointwise ellipticity time diagnostic highlights the same material vortex 
regions labelled as coherent by the other two methods. However, it also suggests further vortical 
regions that are neither rotationally no stretching-wise coherent. The diagnostic does not offer a 
well-defined boundary for the detected vortices either. 

Fig. 16b confirms the expectation that the variationally derived geodesic vortex detection 
method generally labels the same coherent material vortices, but yields smaller vortex boundaries 
due to its uniform stretching requirement. This, coupled with the significantly decreased computa¬ 
tional cost and coding effort, renders LAVD-based Lagrangian vortex identihcation preferable over 
the other two objective methods considered here. 
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